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Abstract 

We describe the method aUowing quantitative interpretation of absorptive images of mixtures of 
BEC and thermal atoms which reduces possible systematic errors associated with evaluation of the 
contribution of each fraction. By using known temperature dependence of the BEC fraction, the 
analysis allows precise calibration of the fitting results. The developed method is verified in two 
different measurements and compares well with theoretical calculations and with measurements 
performed by another group. 
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I. INTRODUCTION 



Analysis of phase transitions offers valuable data of many physical systems. This is 
particularly true for studies of Bose- Einstein (BE) condensation in diluted gases . While 
such studies are very important scientifically, they pose many experimental challenges. The 
main difficulty is caused by extremely low temperatures in which Bose-Einstein condensates 
(BECs) are created and investigated, on the order of 100 nK. It is therefore essential to 
develop reliable detection/imaging methods for ultracold atoms close to the phase-transition 
point (density and temperature). Several such methods have been developed by many groups 
both destructive, hke absorptive imaging, and nondestructive, like phase contrast |3| and 
polarization imaging ^. The simplest and most often applied detection technique is the 
absorptive imaging and this is the method on which we concentrate in this paper. 

In the absorptive imaging one records images of the shadow cast onto a camera sensor 
by atoms usually released from a trap during their free gravitational fall. This yields 2D 
distribution of optical density which reflects spatial density profile of the atomic cloud. Anal- 
ysis of such profiles allows derivation of all relevant physical parameters of the investigated 
sample. The main difficulty in such analysis is associated with the fact that at finite temper- 
atures the BEC fraction is always associated with some fraction of thermal (non-condensed) 
atoms. The thermal fraction plays a very important role in the data analysis as it allows 
determination of the cloud temperature. This fraction decreases with the falling tempera- 
ture of the cloud. Each of the fractions has different density distribution and contributes 
differently to the recorded image. The coexistence and overlap of the two fractions results in 
a bimodal distribution of the optical density which raises interpretational problems. Below, 
we show that a simplistic analysis of such bimodal distributions by fitting them to a sum 
of the Gaussian and Thomas-Fermi functions corresponding to the thermal and condensate 
fractions, respectively, is not satisfactory and leads to systematic errors. 

The problems associated with the anal_y_sis^of absorption profiles are not new and were 



already noticed by several other groups (5|, la, lD]. One attempt which partly avoids the 
problem of the bimodal distribution is the spatial separation of thermal and BEC fractions 
by the Bragg diffraction or by application of an optical lattice [10|]. These methods, 

however, are neither easy nor ideal as they also can introduce inter-fraction interaction. So 
far, the best known method is to analyze the thermal fraction only in its outer regions. 
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outside the degenerate regime, where it is possible to use a simple classical description j^. 
The main disadvantage of this method is that it is not easy and rather arbitrary to find 
the correct size of the excluded central region. As this region is enlarged, the systematic 
error of the fitted parameters is reduced, but the S/N ratio decreases. On the other hand, if 
the excluded central region is too small, the fit is performed to the data which sample also 
the edges of the degenerate region. The resulting systematic errors can be minimized, for 
example by applying corrections based on a numerical solution of the ideal Bose-Einstein 
distribution 6|, which is not a trivial task. 

The present paper introduces the method allowing quantitative analysis of the absorptive 
pictures which ensures correctness of the size of the excluded degenerated region. We present 
the algorithm for analysis of the bimodal distributions which yields accurate ratio of the 
BEG and thermal fractions at finite temperatures. The analysis allows calibration of the 
thermal fraction fits and minimizes number of measurements necessary to obtain statistically 
meaningful averages. 

Section [III presents the procedure of fitting the bimodal distribution of optical density 
and the method for the fit optimization. In Section IIIII we compare our method with other 
commonly used approaches. In Section [IV] we present examples of analysis of the experi- 
mental data with the two methods, the simplest one which uses a sum of a Gaussian and 
Thomas-Fermi distributions and the one we have developed. The paper is concluded in 
Section |Vl 



II. METHOD FOR ANALYZING THE IMAGES OF A CONDENSATE IN NON- 
ZERO TEMPERATURE 

This Section describes the main points of our method for analyzing the BEG pictures 
and its calibration. 



A. Fitting to the bimodal distribution 

Two-dimensional picture of a column optical density (OD) contains information on the 
spatial distribution of the column atomic density in a cloud n (r, z) = OD (r, z) /o"o, where 
r, z are the radial and axial coordinates, respectively, (Fig.[T]) and ctq = 3A^/2'7r is the normal- 
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ized cross-section for atomic absorption at wavelength A. By column densities we understand 
regular densities integrated over the local sample thickness, i.e. the column length. From the 
OD distribution recorded with light intensity /, the non-saturated distribution ODn can be 
calculated as ODn = OD + (1 — exp{—OD)) I/Isat, where Igat is the saturation intensity for 
the imaging transition. If the expansion of the cloud is small, we have to take into account 
that the cloud can be completely dark for the absorption probe beam (see, e.g. 

Well above the critical temperature, the density distribution in the thermal cloud can be 
described with the classical Boltzmann distribution. The column density is described then 
by the Gaussian function: 



ODGauss{r,z) = ODcpeak exp 



1 / r 
2 



Z - Zc 



(1) 



with (Jz being the half- width of the atomic density distribution in the radial and axial 
directions, respectively, ODcpeak denotes the maximum value of the thermal fraction density, 
and [vc, Zc) are spatial coordinates of the maximum. For temperatures close to and lower 
than the critical value, the density distribution becomes predominantly the Bose distribution. 
Then, if the chemical potential is set to zero, the column optical density can be described 
by the, so called, Bose-enhanced Gaussian function 
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(2) 



where g2 (x) = Yj'^=i i.^"^) I (^^) (see, e.g. fl 

With increasing distance from the position of the maximum density, the series terms in 
numerator of ([2]) decrease to zero. At appropriate distance, function ([2]) becomes the Gauss 
function ([1]) which justifies description of the density distribution at the edges of the thermal 
fraction by function ([T]). Nevertheless, more accurate results are obtained if the first three 
terms of series (E]) are used instead. Including of yet more terms only increases computation 
time without noticeably improving the accuracy. 

In the BEG fraction, on the other hand, the distribution of the column optical density of 
atoms in the Thomas-Fermi regime can be described by the TF profile, a clipped parabola. 



ODrF^r, z) = ODrPpeak max 



0, 1 



r — r. 



Rr 




(3) 
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where Rz, Rr are the TF radii in the radial and axial directions, respectively, and ODxFpeak 
denotes maximum of the condensate optical density. When a condensate is not in the TF 
regime, the density distribution is well approximated by a Gauss function j^. 

All pictures in our procedure are taken with the condensates that were expanding for 
time t after their releasing from MT. Knowing spatial distribution of the atomic density in 
the falling cloud, the initial atomic temperature can be determined jisl : 

where = Uit for i = z,r and Tj = {m/2kB) [ujfaf{t)l{l + ojff^)] are the effective temper- 
atures measured after expansion. Moreover, the distribution allows determination of the 
numbers of atoms in each fraction, the thermal fraction Nth and the condensate fraction 
Nbec'- 

Nth = i2n f'^^^a'^{t)aM (5) 
Nbec = ^^i^^^^^RlmM- (6) 

15 (To 

The values of the distribution parameters are to be derived by fitting the OD distribution 
functions to the experimentally recorded profiles. For pictures corresponding to homogenous 
samples consisting exclusively of either thermal atoms or pure BEC, the functions ([1]), ([2]) 
or ([3]) can be fitted, as appropriate. Such fits are performed for the radial and axial sections 
independently. The sections cross the center of mass which coincides with the maximum of 
optical density. In the least-squares fitting procedure the MINUIT library [l^ was used. 

In a bimodal atom cloud, containing both the BEC and thermal fractions, the recorded 
pictures consist of two regions, the external region occupied by the thermal cloud only and 
the internal one where the two fractions coexist. In the later region and close to the border 
between the two regions the density distribution is distorted by the interaction between the 
fractions and by the Bose enhancement of the thermal fraction 5|, LZj . To reduce the effect 
of this distortion, the fitting has to be performed in several steps. 

The first step is to determine the region occupied by the condensate and its direct neigh- 
borhood. For this sake, we approximate the bimodal distribution by the sum of functions 
([1]) and ([3]) with some offset and fit it to the column density picture of the atomic cloud. 
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FIG. 1: Illustration of the spatial extension of the thermal (outer) and BEG (central) fractions 
fitted by different distributions. The scaling factor S is chosen for each picture individually to 
account for the distortion of the density distribution in the intermediate region, as described in 
Subsec. HTBl 

The fit is performed for each direction independently with the least-squares fitting procedure 
based on the MINUIT library. 

Next, using the parameters derived from the fitting curve ([3]), the initial BEC extension, 
i.e. the TF radii, Rz,Rr, are determined. This allows to subtract the BEC contribution 
from the analyzed picture. To account for distortions in the intermediate region, the size 
of the subtracted area is taken with some margin such that the area dimensions are bigger 
than Rz and Rr by scaling factor S (see Fig. [1]). The procedure of exact determination of 
the S value is described in the following subsection. 

After removing the BEC contribution with appropriate safety margins determined by the 
scaling factor, the remaining image consists already of a pure thermal fraction and can be 
fitted by first three terms of series ([2]) with some background. We do it with the least- 
squares method by the 2D NonLinearFit function in Mathematica 5.1. The fit parameters 
allow calculation of the atom number in a thermal cloud and its temperature. 

Having determined the optical density distribution in the thermal fraction and the back- 
ground level, they can be subtracted from the initial picture with full bimodal distribution. 
Additionally, at this stage data points which are below 5% threshold are rejected to eliminate 
the contribution of the intermediate, distorted region of the cloud picture. 
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To such evaluated data the TF profile ([3]) is fit by the 2D least-squares method, which 
yields atom number in the BEC fraction and the appropriate TF radii. For minimization of 
X^, the NMinimize function in Mathematica is used. 

B. Calibration of the thermal fraction region 

The value of coefficient S affects the calculated temperature and atom number of the BEC 
fraction in a bimodal distribution. Taking too big margin, i.e. too big S, eliminates too 
large region which contains information on the density distribution of the thermal fraction, 
thereby lowering the signal/noise ratio and reducing the fit accuracy. Too small value of S, 
on the other hand, introduces systematic errors by including the distorted regions. 

Fig. [2] illustrates the effect of the value of scaling factor 5* on the calculated temperature 
of the bimodal cloud in the radial and axial directions. All fits used for this figure were 
performed for the same picture of a thermal condensate. Fig. [2] illustrates that there exists 
a fairly wide range of the S values where the determined temperature does not change by 
more than one standard deviation (shaded range in Fig. [2]). For S < 1, the determined 
temperature is underestimated by including the intermediate border region affected by the 
BEC fraction, while for S > 1.8, the temperature is not correctly estimated because of low 
S/N ratio of a too much reduced picture. The choice of appropriate value of S requires taking 
into account also the effect of 5* on the number of atoms in the BEC fraction derived from 
the fit. This effect can be seen when studying the phase-transition plot, i.e. the dependence 
of the BEC size normalized to all atoms in a bimodal cloud, Nq/N = Nbec / {^bec + Nth)-, 

1 /3 

on the reduced temperature T/Tc{N) with Tc{N) = h/ks [N ■ {uj'^uJz/1.202)] being the 
critical temperature. 

Nonetheless, the check that the measured temperature of the cloud is independent of 
the size of the exclusion region is not the sufficient criterion of the fit quality. Particularly, 
if the measured data is heavily affected by noise or if the image sizes of either fraction 
are comparable with the image resolution, the temperature stability region (shaded area in 
Fig. [2]) can be very narrow or even vanish completely. We have, therefore, studied further 
consequences of various choices of the S values. 

In Fig. [3] we depict number of atoms in the BEC fraction versus the reduced temperature 
for four typical images of the bimodal cloud. Bimodal distributions corresponding to different 
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FIG. 2: Temperature values in the radial and axial directions determined for experimental recording 
analyzed with different values of S. The marked region indicates the range of 5 values for which 
the determined temperature does not change by more than one standard deviation. 
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FIG. 3: Condensate fraction versus normalized temperature. The points represent results of fitting 
the bimodal distributions obtained with different values of scaling parameter S (between 0.8 and 
1.8) to four different images (marked by different symbols: ^,Oi+) ^)) described in text. 
The elliptical contours depict the regions of concentration points, i.e. those where the derived values 
of atom number and reduced temperature weakly depend on S. The solid line represents function 
Nq/N = 1 — {T /Tc{N)Y [l5] and the broken line represents behavior of a trapped, semi-ideal Bose 



gas 



values of S were fitted to each of the pictures as described above (Subsec. Ill Ap . For each 
recorded image the fitting procedure was performed for different values of S which were 
increasing by constant increments from 0.8 to 1.8. For a specific range of S the derived atom 
numbers, as well as the temperatures, concentrate around specific values, the " concentration 
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FIG. 4: Examples of the heavily fringed and defringed images. The interference fringes have been 
removed from the original image as described in text. 



points". They change very httle with 5* by no more than ±5%. This fact indicates that 5" 
values within the given range provide optimal separation of the perturbed thermal fraction 
form the non-perturbed one which allows proper description of the thermal cloud by function 
([T]) without sacrificing the S/N ratio too much. 

A very convincing verification of the fit quality is the position of the "concentration 



points" on the phase transition plots relative theoretical curves, like in Fig. [3] 15|, Il6 |. 
If a given "concentration point" appears far from the theoretical curve, it indicates that 
the corresponding image was too much affected by some nonstatistical noise, e.g. caused 
by interference fringes or systematics. In such case, another fitting should be tried with 
another (bigger or smaller) background around the atomic cloud. Our experience shows 
that in about 90% of all cases a single repetition of the fitting procedure yields a good 
result. The remaining 10% is most often associated with systematic errors. 

Interference fringes can be removed to a large extent from the image by a sequential 
subjecting the data to FFT, the mask corresponding to the fringe frequencies and to reverse- 



171]. Examples of the heavy fringed 



FFT. This can be done, e.g. with the ImageJ software 
and defringed images are presented in Fig. HI 

Under conditions of our experiment, the optimum S values are between 1.1 and 1.4 and 
depend on the size of the BEG fraction. The center of the "concentration point" on the 
phase transition plots gives the correct value of S, i.e. the correct size of the excluded de- 
generated region, while the size of the "concentration point" allows to estimate its statistical 
uncertainty. 
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FIG. 5: Expanded dependence of the condensate fraction versus normalized temperature. The 
points represent results of fitting the bimodal distributions to the image generated by simulation, 
as described in text, with the S values in the range between 0.8 and 3.0. The elliptical contour 
marks the regions of concentration of nine points, where the derived values of the atom number 
and reduced temperature weakly depend on S. The solid line represents function Nq/N = 1 — 
{T /Tc{N)Y [isl, the broken line represents behavior of a trapped, semi-ideal Bose gas 




III. COMPARISON TO OTHER METHODS 



In order to compare our method with the other approaches we have analyzed a condensate 
image created by a computer simulation. The simulation created a 2D atomic density 
distribution of a bimodal cloud with 0.4 BEC fraction in a given trapping potential. It 
was based on modelling of the BEC part by the Castin and Dum theory 18|], the thermal 
fraction by the Bose-enhanced distribution (2), and by taking into account their ballistic 
expansion within 22 ms. Finally, a noise typical for the recorded absorption pictures was 
added to such constructed simulation. In our simulation the thermal component is treated 
as an ideal Bose gas, while the condensate part is assumed to be in the TF regime. 

In Fig. we depict the number of atoms in the BEC fraction versus the reduced tempera- 
ture obtained by fitting the bimodal distributions to the image generated by the simulation. 
The S values used in the fits are in the range between 0.8 and 3.0 with a step of 0.1. Simi- 
larly as in Fig. [3l the elliptical contour, indicates the regions of concentration of nine points, 
where the derived values of atom number and reduced temperature weakly depend on 5*. 

The basic results of our analysis, i.e. the number of atoms in both the condensate Nq, and 
thermal Nth-, components, temperature T and the Thomas- Fermi radii Rz,Rr, are compared 
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FIG. 6: Condensate fraction versus normalized temperature - the results generated by different 
methods of image analysis applied to the simulated picture of a bimodal distribution. The value 
used in the simulation is depicted as x. Our method gives the value represented by +, a simple 
2D sum of the Gaussian and TF distributions gives Kl, while a Gauss function used instead of the 
Bose-enhanced Gaussian gives □. The solid and broken lines have the same meaning as in Fig. [5j 




with the values preset in the simulation and generated by other methods. The first method, 
due to its simplicity being probably the most common, fits a sum of the 2D Gaussian ([1]) 
and TF distributions to the experimental data. In the second method the Gaussian 
function is fitted to the wings of the spatial distribution, then subtracted from the whole 
distribution and the remaining data is eventually fitted by the TF function. Fitting a 2D 
Gaussian to the wings of the distribution was widely used in the early experiments on BEG 



[e.g. 



0,0, 



20 



2l|). The noise added to the simulated picture reproduces real experimental 



conditions and causes that neither of the methods gives the perfect fit. Still, the described 
method provides the the closest agreement with the simulation parameters. 



IV. TYPICAL EXAMPLES 



In this Section we discuss results of the analysis of the typical experimental images 
obtained with our setup [22;]. The experiment was devoted to studies of the free-fall dynamics 
of a finite-temperature condensate of *^Rb in the F=2 state and is described in more detail 
elsewhere 23|. Here, we present two examples showing how the results derived with the 



method described in Sec. |TT]compare with those obtained with the simple "2D SUM" method 
based on summation of the Gauss and TF distributions. This comparison well illustrates 
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No 


Nth 


T 


Rz 


Rr 


Simulation 


1 


1 


1 


1 


1 


Our method 


1.0104 


1.0660 


1.0304 


1.0045 


0.995 


2D SUM 


1.0199 


1.0082 


0.8737 


1.0081 


0.995 


GTF 


1.0264 


1.0310 


0.9143 


1.0084 


0.996 



TABLE I: Results of the analysis of the simulated image by three different methods: our method, 
"2D SUM" - the method with fitting a 2D sum of the Gaussian ^ and TF ^ distributions, 
" GTF" - the fitting of the Gaussian function to the wings of the distribution and the TF function 
to the condensate part. 

the potential of our method. 



A. Dependence of Nq/N on T/Tc 

Fig. [7] represents a typical experimental dependence of the BEC fraction Nq/N on the 
reduced temperature, T/Tc, analyzed with two approaches. The points represented in Fig. 
[Tl^a) are obtained by using the 2D fitting of the sum of functions ([1]) and ([3]) to the sections 
of absorptive images. The points in Fig. El^b) are obtained using our fitting procedure (Sec. 
HTll . As before, the solid lines are the Nq/N = 1 — {T/Tc{N)Y functions, while the broken 



ones represent results of the calculations along the lines of Ref. [16|. According to Ref. 
jsl, the bigger is the BEC fraction in the sample, the more the results of a simplistic fit 
with a sum of functions ([1]) and ([3]) deviate from the real temperature. The experimental 
data points are obtained from relatively small number of 210 BEC images taken at different 
temperatures. As can be seen in Fig. [7](b), experimental points evaluated with our method 
are in excellent agreement with the model of Ref. while those shown in Fig. [7] (a) deviate 
dramatically from theoretical predictions. 



B. Temperature dependence of the aspect ratio Rr/Rz of a free falling BEC 

Fig. [8] presents results of our measurements of the BEC aspect ratio as a function of the 
reduced temperature, T/Tq, evaluated from 150 images of BEC taken at different temper- 
atures after t = 15 ms free fall. As described previously, the data points in Fig. [8]^a) were 
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TATc T/Tc 



FIG. 7: Comparison of the typical experimental dependences of the BEC fraction, Nq/N, on the 
reduced temperature, T/Tc{N), yielded by 210 images. The points in (a) result from simple fitting 
the 2D sum of functions ([1]) and ([3]) to the sections of absorption images, whereas the points in 
(b) correspond to the fitting procedure described in this paper. The solid and broken lines are the 
same as in Figs. [3l O and[6l 



obtained by 2D fitting of a sum of functions ([T]) and ([3]) to the sections of absorptive images 
whereas those in Fig. MJo) by using our new fitting method. The points evaluated with the 
new method behave qualitatively in the same way as in a similar experiment of Gerbier 
et al. jol- However, in Ref. the thermal and BEC fractions were completely separated 
spatially by Bragg diffraction which eliminated problems of their proper identification in the 
absorptive images. Despite different methods of the fraction separation, our method yields 
qualitatively similar results 2^. On the other hand, the points those obtained with the 
simple "2D-SUM" method exhibit distinctly different, nonphysical behavior . 



V. CONCLUSIONS 



We have developed the method allowing proper interpretation of absorptive images of 
mixtures of BEC and thermal atoms which reduces possible systematic errors arising from 
non-Gaussian distribution of ultra-cold thermal atoms. 

The developed algorithm is based on the fitting procedure of 2D density distributions to 
the absorption profiles describing the thermal fraction. By using the well known tempera- 
ture dependence of the BEC fraction, the analysis allows precise calibration of the fitting 
results and, consequently, reduces number of measurements necessary to obtain statistically 
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FIG. 8: Comparison of the typical dependences of the BEC aspect ratio on the normahzed tem- 
perature, T/Tq, obtained from 150 images taken after t = 15 ms free fall, (a) depicts results of 
2D fitting of a sum of functions ([T|) and ([3]) to the sections of absorption images, points in (b) are 
obtained by application of the described procedure. 



meaningful average values. We compare our method with the others commonly used. We 
have performed experiments verifying the developed method in two different measurements. 
Comparison of the results analyzed with our method and with the simplest fitting proce- 
dure demonstrates that the described method yields far better accuracy and is less prone 
to systematic errors. The results interpreted with our approach are also consistent with 
theoretical calculations and with the results of measurements performed by another group. 
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